function [theta,resUnitFree,Yfit] = OLSforProj(Y,xGrid,appOrder,scaleX_in_g)
nGrid = size(xGrid,1);
nx    = size(xGrid,2);
if appOrder == 1
    % First order approx
    X = [ones(nGrid,1) xGrid./repmat(scaleX_in_g',nGrid,1)];
elseif appOrder == 2
    % Second order approx
    dimX  = (1 + nx + nx*(nx+1)/2);
    X     = zeros(nGrid,dimX);
    X(:,1:1+nx) = [ones(nGrid,1) xGrid./repmat(scaleX_in_g',nGrid,1)];
    idx = 1+nx;
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            if i == j
                X(:,idx) = xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1);
            else
                X(:,idx) = 2*xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1);
            end
        end
    end
elseif appOrder == 3
    % Third order approx
    dimX  = (1 + nx + nx*(nx+1)/2 + nx*(nx+1)*(nx+2)/6);
    X     = zeros(nGrid,dimX);
    X(:,1:1+nx) = [ones(nGrid,1) xGrid./repmat(scaleX_in_g',nGrid,1)];
    idx = 1+nx;
    for i=1:nx
        for j=i:nx
            idx = idx + 1;
            if i == j
                X(:,idx) = xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1);
            else
                X(:,idx) = 2*xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1);
            end
        end
    end
    for i=1:nx
        for j=i:nx
            for k=j:nx
                idx = idx + 1;
                if i == j && j == k
                    X(:,idx) = xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1).*xGrid(:,k)/scaleX_in_g(k,1);
                elseif i == j && j ~= k
                    X(:,idx) = 3*xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1).*xGrid(:,k)/scaleX_in_g(k,1);
                elseif i ~= j && j == k                    
                    X(:,idx) = 3*xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1).*xGrid(:,k)/scaleX_in_g(k,1);
                elseif i ~= j && j ~= k
                    X(:,idx) = 6*xGrid(:,i)/scaleX_in_g(i,1).*xGrid(:,j)/scaleX_in_g(j,1).*xGrid(:,k)/scaleX_in_g(k,1);
                end                
            end
        end
    end
end
Y           = Y./(1+sum(xGrid.^2,2));
X           = X./repmat(1+sum(xGrid.^2,2),1,idx);
theta       = ((X'*X)\(X'*Y))';
Yfit        = (X*theta')';
resUnitFree = Y./Yfit'-1;